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A  polymer-electrolyte  fuel  cell  model  that  incorporates  the  effects  of  hydrogen  sulfide  contaminant  on 
performance  is  developed.  The  model  is  transient,  fully  two-phase  and  non-isothermal  and  includes  a 
complex  kinetic  mechanism  to  describe  the  electrode  reactions.  Comparisons  between  the  simulation 
results  and  data  in  the  literature  demonstrate  that  known  trends  are  well  captured.  The  effects  of  temper¬ 
ature  and  relative  humidity  variations  in  the  anode  stream  are  investigated,  with  further  comparisons  to 
experimental  data  and  a  proposed  explanation  for  the  nonlinear  behaviour  observed  in  the  experiments 
of  Mohtadi  et  al.  [R.  Mohatadi,  W.-K.  Lee,  J.  van  Zee,  Appl.  Catal.  B  56  (2005)  37-42)].  Extensions  to  the 
model  and  future  work  are  discussed. 
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1.  Introduction 

The  proton  exchange  membrane  (PEM)  fuel  cell  has  significant 
potential  as  an  efficient  and  environmentally  friendly  power  source. 
The  main  electrochemical  reactions  involve  hydrogen  and  oxygen 
and  in  the  absence  of  additional  “side”  reactions  the  net  product  is 
simply  water.  However,  the  main  intended  sources  of  hydrogen  (H2) 
are  reformed  hydrocarbons  [2],  which  contain  traces  of  ammonia, 
hydrogen  sulfide  and  carbon  monoxide.  Each  of  these  compounds 
engender  reactions  that  seriously  degrade  performance  and  can 
potentially  cause  long-term  harm  to  the  positive  electrode  [3-5]. 

Carbon-monoxide  (CO)  poisoning  has  been  studied  quite  exten¬ 
sively  and  the  highly  detrimental  effects  of  small  traces  of  CO  in 
the  anode  fuel  stream  are  well-documented  [3,6,7],  Several  mitiga¬ 
tion  techniques  have  been  proposed  and  analysed,  notably  oxygen 
bleeding  and  the  use  of  alternative  catalysts  such  as  platinum- 
ruthenium.  The  influence  of  hydrogen  sulfide  (H2S),  although  first 
studied  as  early  as  1971  by  Loucka  [4],  is  less  well  characterised.  Sev¬ 
eral  researchers  have  shown  that  small  traces  of  H2S  severely  affect 
performance,  on  a  scale  at  least  comparable  with  CO  [1,4,8-18], 
However,  few  mitigation  strategies  have  been  proposed  and  it  has 
been  reported  that  platinum-ruthenium,  while  effective  against  CO, 
is  not  H2S  tolerant  [13,15].  There  is  a  clear  need,  therefore,  for  a 
better  fundamental  understanding  and  characterisation  of  the  H2S 
poisoning  process. 
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Modelling  and  simulation  are  well-established  tools  in  PEM 
fuel  cell  research,  particularly  with  a  focus  on  the  effects  of  water 
retention  in  the  cathode  at  steady  state.  The  methodologies  and 
techniques  were  recently  reviewed  in  [19,20],  Examples  of  time- 
dependent  models  can  be  found  in  [21-23]  and,  in  contrast  to 
steady-state  models,  are  relatively  few  in  number,  even  though 
in  potential  dynamic  applications  such  as  automobile  power,  cell 
stacks  rarely  operate  at  steady  state.  Degradation  mechanisms, 
moreover,  are  inherently  transient,  impacting  performance  on  a 
long  time  scale,  as  in  the  case  of  carbon  corrosion,  or  on  a  relatively 
short  time  scale,  as  with  H2S  poisoning.  Models  of  such  mecha¬ 
nisms,  to  which  one  can  add  radical  attack  of  the  membrane  and 
platinum  dissolution  and  sintering,  are  rare.  This  is  somewhat  sur¬ 
prising  considering  that  the  commercial  viability  of  PEM  fuel  cells 
is  largely  dependent  upon  overcoming  or  minimising  degradation 
phenomena,  many  of  which  can  potentially  be  simulated  on  a  com¬ 
puter  in  a  fraction  of  the  time  required  for  long-life  experiments. 

It  is  also  important  to  note  that  degradation  is  typically  influ¬ 
enced  strongly  by  the  properties  of  the  cell  components,  by  heat 
and  mass  transport,  and  by  the  operating  conditions— this  is  cer¬ 
tainly  true  of  CO  poisoning  and  membrane  failure.  In  order  to  keep 
fitting  parameters  to  a  minimum  and  simulate  operation  over  a 
broad  range  of  conditions  a  model  should  ideally  include  these  fea¬ 
tures  explicitly.  An  earlier  paper  laid  the  foundations  for  a  transient 
model  of  CO  poisoning  and  oxygen  bleeding  on  both  platinum  and 
platinum-ruthenium  that  explicitly  incorporated  the  entire  MEA, 
two-phase  flow  and  temperature  variations  [24],  In  this  paper  a 
similarly  detailed  model  is  developed  to  study  the  H2S  poison¬ 
ing  mechanism,  using  the  results  in  [1,4,8-18]  to  derive  a  detailed 
sub-model  for  the  anode  kinetics. 
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Nomenclature 

A  specific  surface  area  of  agglomerates  (m-1 ) 

a  specific  surface  area  of  platinum  (m-1 ) 

aw  water  activity 

c  molar  concentration  (mol  m-3 ) 

C  specific  heat  capacity  (J  kg-1  K-1 ) 

d  mean  pore  diameter  (m) 

D  diffusion  coefficient  (m* 2 3 4  s-1 ) 

E  cell  voltage  (V) 

Eo  open  circuit  voltage  (V) 

F  Faraday’s  constant  (C  mol-1 ) 

h  mass  transfer  coefficient  (s-1 ) 

H  Henry’s  constant 

i  interaction  parameter  (J  mol-1 ) 

j  current  density  (Am-2) 

J  Leverette  function 

k  thermal  conductivity  (W  m-1  K-1 ) 

L  thickness  (m) 

m  loading  (kg  m-2) 

M  molar  mass  (kg  mol-1 ) 

N  agglomerate  density  (m-3) 

p  liquid  pressure  (Pa) 

pg  gas  pressure  (Pa) 

q  surface  reaction  rate  (mol  m-2  s-1 ) 

tfo2  ORR  rate  (mol  m-3  s-1 ) 

r  reaction  rate  constants  (reaction  dependent) 

R  molar  gas  constant  ( mol  m-3  s-1 ) 

Ragg  agglomerate  radius  (m) 

s  saturation 

S  source/sink  (mol  m-3s-1) 

t  time  (s) 

T  temperature  (K) 

v  velocity  (ms-1) 

x  mole  fraction 

y  distance  (m  or  p,m) 

Greek  letters 

a  charge  transfer  coefficient 

y  diffusion  rate  through  films  (s-1 ) 

5  film  thickness  (m) 

As  entropy  (J  mol-1  K-1 ) 

6  volume  fraction 

q  overpotential  (V) 

6  surface  coverage 

k  absolute  permeability  (m2 ) 

X  membrane  water  content 

pi  dynamic  viscosity  (kg  m-1s-1) 

v  fixed  charge  site  concentration  (mol  m-3 ) 

p  density  (kg  m-3 ) 

a'  surface  tension  (N  m-1 ) 

a  conductivity  (S  m-1 ) 

ct,  collision  diameter  for  species  i  (A) 

r  molar  area  density  (mol  m-2 ) 

0  potential  (V) 

X  contact  angle  (deg) 

S 2  collision  integral  (m-1 ) 

Subscript 
a  anode 

agg  agglomerate 

b  backward  (reaction) 

c  cathode 


C  catalyst  layer 

cap  capillary 

d  dissolved 

e  electrolyte 

g  gas 

G  gas  diffusion  layer 

/  forward  (rate  constants) 

v  vapour 

1  liquid 

p  pore  space 

pt  platinum 

ref  reference 

s  solid/electronic 

0  reference 

v  1  vapour  to  liquid 
v  d  vapour  to  dissolved 
d  -o- 1  dissolved  to  liquid 

Superscript 

l  equilibrium  between  liquid  and  dissolved  phases 

*  equilibrium  between  vapour  and  dissolved  phases 

-  boundary/initial  value 

volume  average 


The  next  section  describes  the  conservation  equations  and 
kinetic  model.  In  Section  3  simulation  results  are  presented,  with 
a  comparison  to  experimental  data  available  in  the  literature.  The 
effects  of  channel  temperature  and  water  activity  variations  on  the 
poisoning  process  are  then  investigated,  with  further  comparisons 
to  experimental  data  from  Mohtadi  et  al.  [1].  A  summary  of  the 
results  and  a  discussion  of  future  work  are  given  in  Section  4. 

2.  Model 

In  this  section  the  main  features  of  the  model  and  the  underly¬ 
ing  assumptions  are  enumerated,  followed  by  a  description  of  the 
conservation  principles,  the  kinetic  model,  the  initial  and  bound¬ 
ary  conditions  and  the  numerical  implementation.  Details  of  the 
parameters  and  the  fitting  procedure  are  also  provided. 

1.  Domain.  The  domain  includes  the  entire  MEA  as  depicted  in 
Fig.  1 .  Each  component  is  modelled  explicitly.  Where  convenient 
the  following  notation  is  used:  GDL  for  the  Gas  Diffusion  Layer, 
CCL  for  the  Cathode  Catalyst  Layer  and  ACL  for  the  Anode  Catalyst 
Layer. 

2.  Catalyst  layers.  The  so-called  “agglomerates  model”  is 
employed,  in  which  the  carbon  support  is  assumed  to  form 
spherical  clusters,  surrounded  by  layers  of  electrolyte  and 
liquid  water.  The  pores  between  agglomerates  are  referred  to 
as  primary  pores,  distinct  from  the  smaller  pores  between  the 
carbon  particles.  Sufficient  contact  between  the  agglomerates 
to  permit  electron  and  proton  migration  is  assumed. 

3.  Reactant  transport  and  transfer.  The  reactants  are  considered 
to  exist  as  species  in  both  the  gas  and  electrolyte  phases.  Devia¬ 
tions  from  Henry’s  law  provide  a  driving  force  for  interfacial  mass 
transfer.  Water  exists  in  three  forms:  dissolved,  vapour  and  liq¬ 
uid.  Mass  transfer  is  driven  by  deviations  from  an  appropriate 
equilibrium. 

4.  Water.  Water  is  considered  to  exist  in  three  forms:  as  a  dissolved 
species,  as  vapour  and  as  liquid.  It  assumed  that  the  net  water 
produced  is  in  liquid  form.  Condensation  and  evaporation  are 
modelled  using  the  approach  in  [25-27],  and  references  therein, 
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dictated  by  the  deviation  of  the  local  thermodynamic  state  from 
equilibrium.  In  a  similar  fashion,  phase  change  between  vapour 
and  dissolved  water  and  between  liquid  and  dissolved  water  are 
introduced  by  considering  deviations  from  appropriate  equlib- 
rium  states. 

5.  Charge.  For  the  conservation  of  charge,  electroneutrality  and  a 
pseudo  steady-state  apply.  The  justification  can  be  found  in  [28], 

6.  Kinetics.  The  kinetics  of  oxygen  reduction  are  assumed  to  follow 
a  global  mechanism,  as  is  usually  the  case  in  modelling  studies. 
The  kinetic  model  at  the  anode  is  derived  from  the  experimental 
results  in  [1,10-12,14], 


Sources  and  sinks  for  the  gas  phase  Eqs.  (1 )— (2) 

Term  ACL  CCL  GDL 

Spi  hpejCHiC,  -  cf )  Vi(Hici  -  cf)  0 

Sv«l  -fWlCXvPg-Psat)  -IW(XvPg-Psat)  ~hv„  i(XvPg-Psat) 

Svwd  fWd(<40  -  cH20)  hv«-d(c|520-C*20)  0 

These  terms  represent,  from  top  to  bottom,  reactant  dissolution  in  electrolyte,  con¬ 
densation/evaporation,  and  vapour/dissolved  water  mass  transfer. 


be  determined  using  the  same  formula.  The  collision  diameters  can 
be  approximated  by 

<tn  2,i=  + 


where  or,  are  the  collision  diameters  for  the  individual  species  i.  The 
values  are  given  in  Table  7  together  with  the  collision  integrals. 

The  gas  velocity  is  given  by  Darcy’s  law  for  flow  through 
a  porous  medium  and  the  absolute  permeability  given  by  the 
Kozeny-Carman  law: 


9Pg 
9y  ’ 


d2  e3 

K  *d-6)2’ 


(4) 


where  k  is  the  absolute  permeability  of  the  gas  diffusion  or  catalyst 
layers,  \i  is  the  dynamic  viscosity  of  the  gas,  d  is  a  mean  pore  diam¬ 
eter  and  K  is  the  Kozeny-Carman  constant.  The  source  terms  are 
given  in  Table  1.  Spi  is  the  rate  of  mass  transfer  between  the  elec¬ 
trolyte  and  gas  phases,  Sv«i  is  the  rate  of  condensation/evaporation, 
and  Sv^d  is  the  rate  of  water  mass  transfer  between  the  electrolyte 
and  gas  phases.  In  Table  1  /ipe  f  are  volumetric  mass-transfer  coef¬ 
ficients  from  the  gas  to  the  membrane  or  electrolyte  phase  on  the 
gas  side  and  H,  are  dimensionless  Henry  constants.  Each  hpei  is 
approximated  based  on  a  local  Sherwood  number  of  2  (for  flow 
past  a  spherical  particle  [29]): 


2.J.  Conservation  principles 

2.1.1.  Reactant  mass  conservation 

Let  s  denote  liquid-water  saturation,  and  Cj  and  c^0  the  pore 
concentrations  of  species  i  =  O2,  H2,  N2,  H2S,  SO3  and  vapour, 
respectively.  The  concentration  of  species  i  dissolved  in  the  elec¬ 
trolyte  and  membrane  is  denoted  cf.  Mass  balance  equations  in  gas 
phase  are  derived  by  taking  into  account  transport  by  diffusion  and 
convection,  and  mass  transfer  to  and  from  the  electrolyte: 

iw-sjq)- A  -vgq)  =  -Spi  (1) 

>-*h20)-|  (^o^-v^o)  =  wvs_d  (2) 


hpe,i 


aShDj 

~d~ 


0(105), 


where  aisthespecific  surface  area  of  the  agglomerates.  The  conden¬ 
sation/evaporation  coefficient  hv^t  is  defined  later.  The  quantity 
v  =  1800  mol  m-3  is  the  fixed-charge  site  concentration  of  the 
membrane, 

The  electrolyte  volume  fraction,  ee,  has  two  components:  one 
from  the  films  that  coat  the  agglomerates  (e{)  and  the  other  from 
the  electrolyte  contained  in  the  agglomerate  interiors  (e[,),  with: 


To  account  for  the  volume  change  due  to  swelling,  which  is 
assumed  to  impact  only  the  film  thickness,  the  following  relation¬ 
ship  is  used: 


where  D,  is  the  free-space  diffusion  coefficient  of  species  i  in  the 
pore  space  and  Dj^0  is  the  free-space  diffusion  coefficient  for  water 
vapour;  T is  temperature  and  pg  is  the  gas  pressure;  e  is  the  poros¬ 
ity,  which  takes  the  value  e  =  ep  in  the  catalyst  layers  and  e  =  €c 
in  the  gas  diffusion  layers.  The  Chapman-Enskog  approximation 
[29],  and  Bruggemann  corrections  have  been  used  for  the  diffusion 
coefficients,  with  values  referred  to  nitrogen: 


D, 


=  0.01858[e(l  —  s)]3/2 


PS  °j2N2^'N2 


(3) 


in  m2  s-1.  For  each  species  i  =  M,  is  the  molar  mass,  erlN2  is  the 
mean  collision  diameter  (Lennard-Jones  force  constant)  and  C2iNl 
is  the  collision  integral.  The  vapour  diffusion  coefficient,  DJ^q,  can 


e{  =  e{0  +  0.0126A, 

where  A  is  the  membrane  water  content  (mol  H2 O/mol  S03-)  and 
e{  0  represents  the  volume  fraction  of  film  without  any  swelling. 
The  latter  quantity  is  related  to  the  film  thickness  without  swelling, 
5e,o. as  shown  in  Table  6  and  derived  in  [23  ].  The  platinum  inside  the 
agglomerates  is  assumed  to  be  inactive  and  therefore  not  to  con¬ 
tribute  to  the  electrochemical  reaction.  Thus  all  reaction  occurs  on 
the  agglomerate  surfaces.  The  combined  volume  fraction  of  the  car¬ 
bon,  platinum  and  small  pores,  ea,  is  assumed  constant.  The  volume 
fraction  of  primary  pores  is  thus  ep  =  1  —  ea  —  ee- 

Mass  balance  equations  for  species  i  =  O2,  H2,  N2,  H2S, 
S03dissolved  in  the  electrolyte  and  membrane  are  derived  by  tak¬ 
ing  into  account  transport  by  diffusion,  mass  transfer  to  and  from 
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Table  2 

Sources  and  sinks  for  the  dissolved  reactants,  potential  and  dissolved-water  Eqs.  (5), 


(13)  ar 

id  (6),  respectively 

ACL 

CCL 

Meaning 

So, 

0 

\qo2 

O2  consumption 

Sh2 

a(-g 3  +  I1J2) 

0 

H2consumption 

Sn2 

0 

0 

n2 

consumption 

Sh2s 

-a(q,  -MJ2+95) 

0 

h2s 

consumption 

Sso3 

aq7 

0 

SO3  production 

Se 

aF[q4  +  2q5  +6q7) 

f<Jo2 

Ss 

-Se 

-Se 

Electron  source 

Sd„l 

hd^ i(cj52o  -  42q) 

o-ch2o) 

Liquid/dissolved 

transfer 

Sw 

“3  ^ 

-1102 

production 
(liquid  phase) 

The  last  term  represents  liquid-water  production  in  Eq.  (8). 


6  =  6g  in  the  gas  diffusion  layers  and  e  =  ep  in  the  catalyst  layers. 
JVfH2o.  Mi  and  p  are  the  molar  mass,  viscosity  and  pressure  of  the 
liquid  water,  respectively,  and  k  is  the  absolute  permeability  of  the 
catalyst  or  gas  diffusion  layers.  The  factor  of  s3  in  Eq.  (9)  is  the 
relative  permeability.  The  source  term  Sw  is  defined  in  Table  2  and 
is  the  rate  of  liquid-water  production. 

By  definition: 

P  =  Pg-Pcap,  (10) 

where  Pcap  is  the  capillary  pressure.  Combining  Eqs.  (8),  (9)  and  ( 10) 
yields: 

m  9s  €KP\  d_  (  3  (  dpcap  ds  _  dpg \  \ 

Mh2o  9t  +  Mimh2o  9y\  \  ds  3y  9y  ) ) 

=  -Svh1  +  VSm+Sw.  (11) 


the  gas  phase,  and  consumption  or  generation  by  reaction: 


d_ 

dt 


(«?) 


9 

3y 


HI) 


=  Si  +  Spi, 


(5) 


where  e  =  €e  in  the  catalyst  layer,  e  =  1  in  the  membrane  and 
are  the  free-space  diffusion  coefficients  in  the  electrolyte  used  with 
Bruggeman  corrections.  The  source  terms  Sj,  i=  02,  Eh,  N2,  H2S, 
S03 ,  are  defined  in  Table  2  and  are  the  rates  of  consumption  or  gen¬ 
eration  of  species  i.  The  surface  reaction  rates  q  and  the  volumetric 
reaction  rate  qo2  will  be  defined  in  Section  2.2. 

The  mass  balance  for  water  dissolved  in  the  electrolyte  and 
membrane,  normalised  with  respect  to  v,  c[^0,  is  derived  by  con¬ 
sidering  movement  by  diffusion  and  drag: 


3_ 

dt 


H 


>)- 


3_ 

3y 


<2o 

0  3y 


,3/2. 

44Fv 


3te\ 

9y  J 


For  a  discussion  on  the  form  of  the  permeability  and  cap¬ 
illary  pressure  the  reader  is  referred  to  [26],  In  this  paper  the 
widely  used  Leverette  function  is  adopted  (found  in,  for  example 
[31]): 


Pcap  =  O-'  COS  XGi/l^l-S), 


for  the  gas  diffusion  and  catalyst  layers,  respectively.  In  these 
expressions  &  is  the  surface  tension  and  )  is  the  Leverette  func- 


J(!)  =  1 .417|  -  2.12|2  +  1 ,262|3 


=  -Sv^d-Sd4>l,  (6) 

in  which  e  =  1  for  the  membrane  and  e  =  6e  for  the  catalyst  layer, 
and  is  the  diffusion  coefficient  of  water  in  the  electrolyte 
subject  to  a  Bruggeman  correction. 

The  following  form  of  D^0  for  Nation  is  taken  from  [30]: 

dh2o  =  u iM*  +  161  e_A-)exp  . 


where  u  1  =  4.17  x  10-8  m2  s-1.  The  source  term  Sv^d  was  previ¬ 
ously  defined  and  Sd4>],  given  in  Table  2,  is  the  rate  of  water  mass 
transfer  between  the  liquid  and  dissolved  phases.  Both  are  dis¬ 
cussed  below.  Note  that  the  water  content  and  concentration  are 
related  by 


‘  u2  -  0.0126c£2O  ’ 


(7) 


where  U2  =  1  molm-3.  Eq.  (7)  will  be  used  in  the  sequel. 

A  mass  balance  of  liquid  water  is  derived  considering  transport 
by  convection,  driven  by  pressure  gradients: 


gpi  9s  9  (  epi  \ 
Mh2o  dt  3y  yMH2o  / 


=  — .S'V( . |  +  vSd^i  -(-  Sw , 


(8) 


where  the  liquid-water  interstitial  velocity  V]  is  given  by  Darcy’s 
law: 


Xc  and  kc  are  the  contact  angle  and  absolute  permeability  of 
the  catalyst  layers,  respectively,  and  xg  and  kq  are  the  con¬ 
tact  angle  and  absolute  permeability  of  the  gas  diffusion  layers, 
respectively. 


2.1.2.  Charge  conservation 

Equations  for  the  potentials  in  the  electrolyte/membrane  and 
carbon,  <pe  and  <ps  respectively,  are  derived  from  conservation  of 
charge  in  each  phase,  assuming  electroneutrality  and  steady-state 
conditions,  as  justified  in  the  assumptions: 


where  ae  and  ers  are  the  protonic  and  electronic  conductivity 
respectively,  used  with  Bruggeman  corrections,  and  F  is  Faraday’s 
constant,  e  =  ee  in  the  catalyst  layer  and  e  =  1  in  the  membrane  for 
<j)e,  while  c  =  1  —  ee  in  the  gas  diffusion  layer  and  6  =  es  in  the  cat¬ 
alyst  layer  for  <ps-  The  source  terms  Ss  and  Se  are  defined  in  Table  2. 
The  protonic  conductivity  is  assumed  to  take  the  form  given  in  [32] 
for  Nafion: 


ffe  =  U3  exp 


/ 1286 
V  303 


1^)(0.514A.-0.326), 


(14) 


Vi  =  — 


ks3  9p 
Mi  9y' 


where  113  =  1  S  m_1.  The  electronic  conductivity  is  constant  and  its 
value  is  given  in  Table  10. 
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Term 

Membrane 

ACL 

CCL 

GDL 

Qact 

0 

aFrj3q4 

FMo2 

0 

Qrev 

0 

-A  saTaq4 

Ascr<jo2 

0 

Qohm 

«*m2 

Qpc 

0 

-hgiVa 

-hgiSv^i 

-hglSv* 

They  a 

re,  from  top  to  bottom,  activation  losses 

i,  heats  of  reaction,  ohmic  losses  and  he, 

ats  of  evaporation. 

2.1.3. 

Energy  conservation 

or,  equivalently: 

A  thermal  energy  balance  equation  is  derived  by  taking  into 
account  heat  generation  and  heat  transport  through  the  gas,  liquid 
and  solid  phases,  and  assuming  a  single  temperature: 


I  espiC|V|T  +  e(l  -s)pgCgvgT 


in  the  catalyst  and  gas  diffusion  layers,  and 


3 1 


(PmCmT)  + 


=YlQk’ 

(15) 

(16) 


in  the  membrane.  In  these  equations  p\,  pg,  pm  and  ps  are  the  den¬ 
sities  of  the  liquid,  gas,  membrane  and  solid  phases,  respectively, 
and  Q,  Cg,  Cm  and  Cs  are  the  specific  heat  capacities  of  the  liquid, 
gas,  membrane  and  solid  phases,  respectively,  k  and  pCp  are  the 
volume-averaged  thermal  conductivity  and  thermal  capacitance: 


k  =  kp{  1  -  s)e  +  k]Se  +  ks(l  -  e), 

pCp  =  espiQ  +  e(l  -  s)pgCg  +  psCs(l  -  e), 


(17) 


where  kp,  ks  and  Jq  are  the  thermal  conductivities  of  the  pore  space, 
solid  (averaged  over  carbon,  electrolyte  and  platinum)  and  liquid 
water,  respectively.  In  the  gas  diffusion  layers  e  =  Cq  and  in  the 
catalyst  layers  e  =  ep.  In  the  membrane  the  only  form  of  heat  trans¬ 
port  is  conduction.  The  gas  phase  thermal  conductivity  and  specific 
heat  capacity  are  approximated  by  values  for  air.  The  heat  genera¬ 
tion  terms  Q*  are  defined  in  Table  3.  In  these  expressions  hgl  is  the 
liquid-gas  enthalpy  change  for  water,  —  Asc  is  the  entropy  associ¬ 
ated  with  the  oxygen  reduction  reaction  and  — Asa  is  the  entropy 
associated  with  the  hydrogen  oxidation  reaction.  Note  that  as  a  sim¬ 
plification,  the  heats  of  reaction  of  the  other  anode  reactions  are 
neglected. 

2.1.4.  Water  phase  change 

The  treatment  of  water  mass  transfer  between  the  three  phases 
is  now  detailed.  Condensation  and  evaporation  are  driven  by  the 
deviation  from  equilibrium:  xvpg  -psa t,  where  psa t  is  the  satura¬ 
tion  pressure  of  water  and  the  first  term,  in  which  xv  is  the  vapour 
mole  fraction,  is  the  partial  pressure  of  the  vapour.  The  condensa¬ 
tion/evaporation  coefficient  hV4¥\  in  Table  1  takes  the  form: 

%(1  -s)xv 

(18) 


!iv. 


hcond  PT  Xvpg-psat>  0 
XvPg-Psat  <  0 


ePSA 

,pMH2o 


(21) 


cft2o(u2  +  °.0126^)  =  **.  (20) 

In  these  formulae,  aw  =  xvpg/psa t  is  the  water  vapour  activity. 
The  mass-transfer  coefficient  hv„d  is  approximated  from  the  results 
in  [34]: 

f  3ads,v(  i  -  s)X  c{520  -  cfj20  <  0 
=  <  . 

[  3desv(l  -  s)X  cH20-cJ20>  0, 

where  hdes  v  and  /tadsv  are  desorption  and  adsorption  coefficients, 
respectively.  Their  values  are  given  in  Table  10. 

The  equilibrium  membrane  water  content  depends  on  its  envi¬ 
ronment,  with  either  relationship  (19)  for  contact  with  vapour  or 
X  =  Xl  =  16.8  for  contact  with  liquid  water.  Note  that  the  liquid- 
equilibrated  dissolved  water  concentration,  c^,  is  given  by  Eq. 

(7): 

j 


H2°  u2  +  0.0126A.' ' 

The  discontinuity  between  the  vapour-saturated  and  liquid  val¬ 
ues  is  known  as  Schroeder’s  paradox.  The  mass-transfer  term  Sd4>1 
in  Eqs.  (6)  and  (8),  and  Table  2,  is  decomposed  into  terms  for  absorp¬ 
tion  and  desorption  of  liquid  water  to  and  from  the  electrolyte  in 
the  catalyst  layer.  When  the  liquid-equilibrated  water  content  value 
is  reached  or  exceeded,  it  is  assumed  that  desorption  of  water 
from  the  electrolyte  takes  place  (as  liquid),  the  magnitude  of  which 
is  driven  by  c^0  —  c[ ^Q.  Adsorption  is  assumed  to  take  place  for 
ch  o  <  ch  o  Provided  s  >  s„,  where  s,  is  the  immobile  saturation. 
The  coefficient  hd^  in  Table  2  therefore  takes  the  form: 

hd«*l  =  W^o  -  c^o)  +  hads,,H(s  -  stm-cii0  +  c'H20),  (22) 

where  H(  )  is  the  Heaviside  function  and  hdesI  and  /iadsd  are  the 
coefficients  of  desorption  and  absorption,  which,  for  simplicity,  are 
assumed  to  be  constant.  Their  values  are  given  in  Table  10. 

2.2.  Reaction  kinetics 

The  reaction  rates  appearing  in  Tables  2  and  3  are  yet  to  be 
defined.  The  kinetics  at  the  cathode  are  approximated  by  the  fol¬ 
lowing  global,  one-step  reaction: 

Oxygen  reduction  reaction (ORR) :  02  +  4H+  +  4e_  ^  2H20, 


/icond  and  hemp  are  the  condensation  and  evaporation  rate  con¬ 
stants,  whose  values  are  taken  from  [25]. 

In  a  similar  fashion,  the  vapour-dissolved  phase-change  term  in 
Eqs.  (2)  and  (6),  Svt^d,  is  driven  by  the  deviation  from  equilibrium 
between  the  vapour  and  dissolved  water,  cf{  Q  -  cf,  0,  where  c*  Q  is 
the  dissolved  water  concentration  at  equilibrium,  as  given  in  [33]: 


9o2(>]c,T,cs02)  = 


e  the  Butler-Volmer  law  (in  mol  m-3  s_1 ): 

(-“#)}■ 


=  0.3  +  10.8aw  -  16 al  +  14.1a^, 


(19) 


The  various  terms  in  the  formula  are  the  exchange  current  den¬ 
sity  jo2,ref:  the  anodic  and  cathodic  transfer  coefficients  a3  and  ac, 
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respectively;  the  reference  oxygen  molar  concentration  Co2,ref;  the 
volumetric  specific  surface  area  of  catalyst  (per  unit  volume  of  cat¬ 
alyst  layer)  a;  and  the  overpotential  r)c.  The  expression 

_  OptUlpt 

a~  I 

relates  the  specific  surface  area  of  platinum  to  the  mass  specific 
platinum  surface  area  (platinum  surface  area  per  unit  mass  of 
platinum),  apt,  the  platinum  loading,  jnpt,  and  the  catalyst  layer 
thickness,  L.  The  overpotentials  in  the  cathode  and  anode,  r)c  and 
ij3  respectively,  are  defined  as 


H2S  +  Pt-S^Pt-S2  +  2H+  +  2e-.  (R6) 

According  to  the  cyclic-voltammetry  measurements  in  [10,13], 
the  main  species  to  form  is  in  fact  Pt-S. 

The  adsorbed  sulfur  can  be  oxidised  via  reaction  with  water, 
as  shown  in  reactions  (R7)  and  (R8)  below.  This  mechanism  was 
proposed  by  Loucka  [4],  and  confirmed  by  Mohtadi  et  al.  [13]  and 
Wang  et  al.  [18],  in  which  the  authors  report  the  presence  of  sulfate 
ions: 

Pt-S  +  3  H20  ^  S03  +  6H+  +  6e_  +  Pt.  (R7) 


3?c  =  0s-0e-Eo,  3?a  =  0s-0e,  (25) 

where  E0  is  the  open  circuit  potential  vs.  S.H.E.  [35]: 

£0  =  1  -23  -  9  x  10“4  (T  -  298.15). 

An  agglomerate  model  is  used  for  the  cathode  catalyst  layer,  in 
which  c*2  is  the  oxygen  concentration  at  the  agglomerate  surfaces. 
It  is  related  to  the  bulk  value,  Co2,  by  balancing  the  rate  of  reac¬ 
tion  with  the  rate  of  diffusion  of  reactant  through  electrolyte  and 
water  films  to  the  surfaces  of  the  agglomerates.  The  final  form  of 
the  reaction  rate  is 

e<*,Fr,c/RT  _  ^-ctcFrjc/RT 

q°2  =  4y(1'C0 2  y+a'(e°‘^c/RT_e-acFnc/RT)  ’  (26) 

where  the  parameter  a'  is  given  by 

,  q/02,ref 

“  “  4fCo2,ref 

y  is  a  measure  of  the  diffusion  rate  through  the  films  [y/a  is  the 
piston  velocity)  and  takes  the  form: 


Pt-S  +  4  H20  ^  SO4-  +  8H+  +  6e“  +  Pt.  (R8) 

The  mechanism  considered  by  Shi  et  al.  [16]  is  an  application  of 
the  general  contamination  model  developed  in  [36],  After  approxi¬ 
mation,  it  consists  essentially  of  the  forward  parts  of  reactions  (R1 ), 
(R3)  and  (R8),  the  competitive  adsorption  of  H2  on  platinum  and 
the  following  reaction: 

Pt-H  +  H2S  -►  Pt-S  +  H2  +  H+  +  e- 

This  reaction,  though  not  supported  by  experiment,  appears  to 
have  been  included  for  consistency  with  the  model  in  [36], 

The  model  developed  below  includes  reactions  (R1  )-(R7).  Reac¬ 
tion  (R6)  is  neglected  as  per  the  results  of  [10,13]  and  reaction  (R8)  is 
neglected  for  simplicity.  It  is  not  expected  that  the  overall  behaviour 
of  the  system  will  be  altered  qualitatively  by  these  assumptions.  The 
model  incorporates  the  two  main  mechanisms  of  adsorption  and 
oxidation,  for  both  hydrogen  and  sulfur. 

Equations  for  the  evolutions  of  the  site  coverages  for  atomic 
hydrogen  and  sulfur,  9h  and  9$  respectively,  are  then  as  follows: 


(A'D^)  (AD«02/8e) 

A'AA+AD^/Se  ’ 

in  which  A'  =  4jr(Ragg  +  Se)2N.  In  these  expressions:  <$e  and  <5|  are 
the  electrolyte  and  liquid-water  film  thicknesses  respectively:  N 
is  the  number  of  agglomerates  per  unit  volume;  A  =  4jzRj&&N  is 
the  specific  surface  area  of  agglomerates,  assuming  that  the  entire 
surface  area  of  each  agglomerate  is  covered;  Ragg  is  the  agglomer¬ 
ate  radius;  and  D|  is  the  diffusion  coefficient  of  02  through  liquid 
water.  The  film  thicknesses,  Se  and  5],  are  defined  in  Table  6  and 
take  account  of  electrolyte  swelling.  The  reader  is  referred  to  [23] 
for  derivations. 

The  H2S  kinetics  on  platinum  were  investigated  by  Contractor 
and  Lai  [10,11],  who  suggested  that  there  are  two  adsorbed  forms 
of  sulfur,  one  strongly  and  one  weakly  bonded.  Mathieu  and  Primet 
[12],  proposed  the  following  form  for  the  surface  reactions: 

H2S  +  Pt^Pt-S  +  H2.  (Rl) 

Pt-H +  H2S^  Pt-S  +  |h2.  (R2) 

These  reactions  occur  in  competition  with  the  adsorption  and 
electro-oxidation  of  hydrogen  [24]: 

H2  +  2Pt  ^  2 Pt-H.  (R3) 


^=-<j2  +  2q3-<j4  (28) 

=qi+q2  +  Qs-Q7  (29) 

where  the  rates  are  defined  by 

41  =  rVcH2s0Pt ^P  ~  ri/,<Vs exp  ^(1 

q2=r2f0nciiS-r2b6s(c^  ' 
q-i  =  rifCdH__elt-r-ih0l 
q4  =  r40H  sinh  ( 

V  RT  )  (30) 

45  =  r5/<4s0Pt  exp  (  -  r5b@SCH+ 

xexp 

47=r7/(^2o)30sexp(^^) 

-r7hcso3CH+0Pt  exp  (~6(1  ~R“7)F??a)  • 

The  quantity  r  is  the  molar  area  density  of  catalyst  sites.  The 
surface  coverage  of  free  platinum  sites,  0Pt,  is  given  by 


Pt-H^H+  +  e-  +  Pt.  (R4) 

Mohtadi  et  al.  [1  ],  suggested  that  platinum  sulfide  (Pt-S)  can  also 
be  formed  electrochemically  according  to  reaction  (R5)  below.  The 
same  reaction  was  proposed  by  and  Najdeker  and  Bishop  [14],  who 
further  suggested  that  platinum  disulfide  (Pt-S2)  could  be  formed 
according  to  reaction  (R6): 

H2S  +  Pt  ^  Pt-S  +  2H+  +  2e“ .  (R5) 


9Pt  =  -i-eH-es  (3i) 

Eqs.  (30)  and  (31 )  are  derived  on  the  assumption  that  coverage 
does  not  exceed  a  monolayer.  The  forms  of  the  reaction  rates  are 
based  on: 

(1)  a  Frumkin  isotherm  for  reaction  (Rl),  sulfur  adsorption  and 
desorption,  with  rate  qi ; 
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(2)  Langmuir  isotherms  for  reaction  (R2)  and  reaction  (R3),  H2 
dissociative  adsorption  and  desorption,  with  rates  q2  and  q3 
respectively; 

(3)  Butler-Volmer  kinetics  for  reaction  (R4),  electro-oxidation  of 
hydrogen,  with  rate  ^4; 

(4)  and  Butler-Volmer  type  kinetics  for  the  oxidation  reactions 
(R5)  and  (R7),  with  rates  qs  and  q 7,  respectively. 

In  the  Frumkin  model  [37],  the  apparent  standard  free  energy  of 
adsorption  is  assumed  to  vary  linearly  with  6S,  with  slope  equal  to 
the  interaction  parameter  i,  [38],  Since  their  values  are  not  known, 
all  charge  transfer  coefficients,  are  assumed  to  be  1/2.  This  is  a 
common  modelling  assumption.  It  was  further  assumed  that  the 
backward  rates  for  reactions  (R5)  and  (R7)  are  negligibly  small. 
The  values  of  ty,  rlb  and  r7f  were  then  fitted  to  achieve  a  qualita¬ 
tive  match  to  experimental  results  in  the  literature.  Similar  fitting 
procedures  have  been  performed  in  [3,6,7,16]. 

2.3.  Initial  and  boundary  conditions 


Table  4 

Operating  conditions  assumed  in  the  calculations,  unless  otherwise  specified 
Symbol  Quantity  Size 


Tc  Cathode  channel  temperature 

Ta  Anode  channel  temperature 

aw,c  Cathode  channel  water  activity 

aw,a  Anode  channel  water  activity 

*o2,c  Oxygen  mole  fraction  in  cathode  channel3 
XH2,a  Hydrogen  mole  fraction  in  anode  channel3 
Xso3.a  S03mole  fraction  in  anode  channel3 

Xn2,c  Nitrogen  mole  fraction  in  cathode  channel3 
XN2,a  Nitrogen  mole  fraction  in  anode  channel3 
CH2s,a  H2S  mole  fraction  in  anode  channel3 

Pa  Gas  pressure  in  the  anode  channel 

Ch2o,c  Vapour  concentration  in  cathode  channel 

CH2o,a  Vapour  concentration  in  the  anode  channel 

Liquid-water  removal  constants  (anode  and  cathode)3 
j'a  Applied  current  density 

E  Cell  voltage 

3  After  subtracting  the  vapour  concentration. 


60  °C 
60  °C 
0.9 
0.9 
0.21 


0.79 

0.6 

2  ppm 
300 kPa 
300 kPa 
6.38  mol  m~3 
6.38  mol  m~3 
0.075  m-1 
0.5  A  cm-2 
0.5  V 


At  the  interfaces  between  the  membrane  and  catalyst  layers, 
y  =  y2  and  y  =  3/3  shown  in  Fig.  1 ,  the  gas-phase  and  liquid-water 
fluxes  are  taken  to  be  zero.  Similarly,  the  fluxes  of  protons  and 
dissolved  species  at  the  interfaces  between  the  catalyst  and  gas 
diffusion  layers,  y  =  y\  and  y  =  y4,  are  negligibly  small: 


y=y\ .  y4 ; 


y  = 


y2,y3  •• 


I,.,  ,3cd 

e3/2Dd_J_  =0 

1  3y 

e3/2ffe^=0 

3y 

nd  9CH20  ,  SAffe  3 0e  _  n 

^20  ty  +  44Fv  3y 


(32) 


(33) 


At  the  interfaces  between  the  channels  and  gas  diffusion  layers, 
the  gas  mole  fractions  are  prescribed  or  calculated  from  the  other 
conditions: 


Xi={kil  y  =  S  ’  '  =  h2,02,H2S,S03,N2,H20.  (34) 


Likewise,  temperature,  water  activity  and  pressure  are  pre¬ 
scribed  according  to  the  channel  values: 


T(y0)  =  Tc  T(y5)  =  Ta, 

Pg(yo)  =  Pc  Pg(ys)  =  Pa,  (35) 

Qw(yo)  =  aw,c  aw(y5)  =  aw,a. 


The  concentrations  of  water  vapour  in  the  cathode  and  anode 
channels,  cv,c  and  cv,a  respectively,  are  calculated  from  the  water 
activities  and  saturation  pressures  ([32]): 

log10Psat  =  -2.1794  +0.02953(7  -  273.15) 

-9.1837  x  10-S(T  -  273.15)2  (36) 

+1.4454  x  10-7(T-273.15)3, 


where  psat,c  is  the  cathode-channel  saturation  pressure  (expressed 
in  Pa).  A  similar  calculation  applies  on  the  anode  side  where  the 
channel  saturation  pressure  ispSat,a- 


For  operation  in  the  potentiostatic  mode  the  cell  voltage,  E,  is 
prescribed  at  the  cathode  channel/gas  diffusion  layer  interface,  and 
at  the  anode  channel/gas  diffusion  layer  interface  the  electronic 
potential  is  assigned  a  value  of  zero: 


0s(yo)  =  £,  0s(y5)  =  O. 


(37) 


When  the  cell  is  operated  in  galvanostatic  mode,  a  current  den¬ 
sity  ja  is  prescribed  at  the  electrode/gas  channel  interfaces: 


a  90s  =  f  ja  y  =  0 
’  p  3y  \  -Ja  y=y5 


(38) 


where  erp  is  the  electronic  conductivity  of  the  current  collector 
(plate). 

The  final  boundary  conditions  are  those  for  liquid  water  at  the 
interfaces  between  the  gas  channels  and  the  gas  diffusion  layer. 


Table  5 

The  default  parameter  values  relating  to  structural  properties 


Symbol 

Quantity 

Size 

I 

Catalyst  layer  thickness 

25  |xm 

Lm 

Membrane  thickness 

tc 

GDL  thickness 

200  Jim 

Electrolyte  volume  fraction  in 
agglomerates3 

0.15 

Volume  fraction  of  agglomerates 
[40] 

0.4215 

Volume  fraction  of  carbon  in  CLb 

0.2 

<+ 

Porosity  of  the  GDL  [41  ] 

0.74 

Ragg 

Agglomerate  radius  [40] 

0.5  jjum 

4,0 

Electrolyte  film  thickness  without 
swelling  [23] 

0.1  jxm 

<4 

Electrolyte  volume  fraction 
without  swelling  [23] 

^[(Ragg  +  4,o)3-Raggl 

4 

Electrolyte  film  thickness  [m]  [23] 

\/Ragg  +  45R  “  Ragg 

8, 

Water  film  thickness  [m]: 

Ri=Ragg  +  4[23] 

\As  +  IHi  _  Ra 

N 

Agglomerate  density11 

5.8  X  1017  m-3 

apt 

Specific  surface  area  of  platinum 
[42] 

1000  cm2  (mg  Pt)_1 

mpt 

Platinum  loading3 

0.4  (mg  Pt)  cm-2 

Xc 

Catalyst-layer  contact  angle  [43] 

90° 

Xc 

Gas  diffusion  layer  contact  angle 
[43] 

120° 

dc 

GDL  pore  size3 

10  (xm 

dc 

Catalyst-layer  pore  size3 

2  M-m 

3  Assumed  value. 
b  Estimated  value. 
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They  are  approximated  using  the  following  steady-state  flux  con¬ 
ditions  aty  =  y0  andy  =ys  (see  [23]  for  details): 

-  PjS  =  0,  j  =  a,c,  (39) 

where  fy  =  0  corresponds  to  zero  water  removal  from  the  anode 
channel,  j  =  a,  or  cathode  channel,  j  =  c. 

The  initial  conditions  for  the  pressures,  temperatures  and 
vapour  concentrations  are  consistent  with  the  conditions  in  the 
channels.  The  electronic  potential  at  the  initial  time  is  given  by  the 
open-circuit  potential  at  the  cathode  and  is  zero  at  the  anode.  The 
initial  protonic  potential  is  zero  in  all  regions.  The  cell  voltage  or 
current  is  ramped  at  a  finite  rate  as  discussed  later.  The  initial  water 
content  of  the  membrane/electrolyte  is  given  by  equilibrium  with 
the  vapour  in  the  channels.  The  initial  liquid-water  saturation  and 
all  surface  coverages  are  uniformly  equal  to  zero  at  the  initial  time. 

2.4.  Numerical  details  and  parameters 

The  initial-boundary  value  problem  developed  above  was  solved 
in  the  software  package  COMSOL  Multiphysics®  on  a  uniform  grid 
(typically  128  points)  using  quartic  Lagrange  polynomials  as  trial 
and  test  functions.  The  relative  tolerance  was  set  to  a  value  of 
1  x  10  6  and  the  absolute  tolerance  to  1  x  10-8.  The  switch  func¬ 
tions  were  substituted  with  hyperbolic  tangent  functions  to  smooth 
the  discontinuities. 

The  default  set  of  parameter  values  is  given  in  Tables  4-10  .  Sev¬ 
eral  parameters  are  estimated,  as  indicated  in  Tables  4-10,  and  the 
rest  are  found  from  the  literature  with  references  provided.  Where 


Table  6 

The  default  parameter  values  relating  to  electrochemical  properties 


Symbol 

Quantity 

Size 

J02,ref 

Cathode  exchange  current 
density  [22] 

10-2  A  m-2 

Reference  02  concentration3 

0.05  mol  m~3 

“c 

Cathodic  charge  transfer 
coefficient 

0.55 

«a 

Anodic  charge  transfer  coefficient 

0.45 

«t 

Charge  transfer  coefficient  for 
reaction  (Rl)3 

0.5 

a4 

Charge  transfer  coefficient  for 
reaction  (R4)3 

0.5 

“5 

Charge  transfer  coefficient  for 
reaction  (R5)3 

0.5 

on 

Charge  transfer  coefficient  for 
reaction  (R7)3 

0.5 

Interaction  parameter3 

10  kj  mol-1 

T 

Molar  area  density  of  platinum 
sites  [7] 

0.01042  mol  m-2 

r'f 

Forward  rate  constant:  reaction 
(Rl)b 

3  x  10-2  ms-1 

rlb/rlf 

Backward  rate  constant:  reaction 
(Rl)3 

1.4  x  10“6 

r2f 

Forward  rate  constant:  reaction 
(R2)b 

1  x  10"2  ms-1 

r2b/r2f 

Backward  rate  constant:  reaction 
(R2)3 

1  x  10-6  m3/2  mol-1/2 

r3 f 

Forward  rate  constant:  reaction 
(R3) [44] 

3e-i0400/i?r  m  s-i 

r3b/r3f 

Backward  rate  constant:  reaction 
(R3) [44] 

4.18  x  1011  e-87900/^  mol  m-3 

r4 f 

Forward  rate  constant:  reaction 
(R4) [44] 

23.1e-1670°/RT  mol  m~2  s"1 

r5 f 

Forward  rate  constant:  reaction 
(R5)3 

1  x  10-4ms_1 

r7f 

Forward  rate  constant:  reaction 
(R7)b 

4  x  10-12  m7  s_1  mol-2 

3  Assumed  value. 


Table  7 

Parameters  for  the  gas-phase  diffusion  coefficient  in  Eq.  (3) 


Symbol 

Quantity 

Size 

ao2 

02 collision  diameter  [29] 

3.433  A 

<*h2 

H2collision  diameter  [29] 

2.915  A 

CTvapour 

Vapour  collision  diameter  [29] 

2.903  A 

ctn2 

N2  collision  diameter  [29] 

3.667  A 

oh2s 

H2S  collision  diameter  [29] 

3.748  A 

ctS03 

S03  collision  diameter  [29] 

4.29  A 

^n2,o2 

02-N2  collision  integral  [29] 

0.966 

^n2,h2 

H2-N2  collision  integral  [29] 

0.848 

lr  Vapour-N2  collision  integral  [29] 

1.305 

flN2,N2 

N2-N2  collision  integral  [29] 

0.949 

fiH2,N2S 

H2S-N2  collision  integral  [29] 

1.129 

£?n2,so3 

S03-N2  collision  integral  [29] 

1.06 

Mo2 

Molar  mass  of  02 

0.0032  kg  mol"1 

MH2 

Molar  mass  ofH2 

0.0002  kg  mol-1 

Mh2o 

Molar  mass  ofH20 

0.0018  kg  mol"1 

Mn2 

Molar  mass  ofN2 

0.0028  kg  mol"1 

Mh2s 

Molar  mass  ofH2S 

0.0034  kg  mol-1 

mso3 

Molar  mass  of  S03 

0.008  kg  mol"1 

Table  8 

Default  parameter  values  related  to  mass  transport 

Symbol 

Quantity 

Size 

°o2 

02  diffusion  coefficient  in  the 
electrolyte1  [45] 

3.1  xl0-7e-2768/Tm2s-) 

°h2 

H2  diffusion  coefficient  in  the 
electrolyte3’ [46] 

6.92  x  10-9  m2  s"1 

°h2s 

H2S  diffusion  coefficient  in  the 
electrolyte3 

4.38  x  10-9  m2  s-1 

topei 

Mass  transfer  ratesb 

105  s-1 

KC  ' 

Absolute  permeability  of  CCL  [47] 

10-13  m2 

Kq 

Absolute  permeability  of  GDL  [47] 

8.7  x  10-12  m2 

Ml 

Liquid  water  viscosity 

10-3  Pa  s 

M 

Dynamic  viscosity:  airb 

1.8  x  10-5  Pas 

M 

Dynamic  viscosity:  H2b 

8.4  x  10_6Pas 

Surface  tension  [26] 

0.07  Nm-1 

<Ts 

Electronic  conductivity 

500  Sm-1 

3  Approximated  by  value  in/for  liquid  water  at  60  °C. 
b  Estimated. 


available,  values  for  Nafion  were  used.  For  several  mass  transport 
and  transfer  parameters,  values  corresponding  to  liquid  water  have 
been  used  as  estimates.  The  pore  size  and  porosity  values  fall  within 
typical  ranges  for  conventional  PEMFC  [39]  and  variations  would 
not  qualitative  affect  the  results  in  this  work. 

The  kinetic  parameters  were  taken  from  the  literature  where 
available.  In  the  absence  of  primary  experimental  data,  the 


Table  9 

Default  parameter  values  related  to  mass  transfer 

Symbol  Quantity 

Ho2  O2  Henry’s  law  constant  [48  ] 

Hh2  H2  Henry’s  law  constant  [49] 

Hh2s  H2S  Henry’s  law  constant 

o3  SO3  Henry’s  law  constant 

hdes,i  Desorption  coefficient  of  dissolved  to 

liquid  water3 

hads,i  Adsorption  coefficient  of  liquid  to 

dissolved  water3 

toads, v  Absorption  coefficient  of  vapour  to 

dissolved  water  [34] 

hdes  v  Desorption  coefficient  of  dissolved 

water  to  vapour  [34] 

toevap  Evaporation  coefficient  [25] 

tocond  Condensation  coefficient  [25] 

X*  Liquid-equilibrated  water  content  [25] 


Size 


0.15 

0.63 

0.32 

1.94 

100 

10 


100  s-1 
100  r1 
16.8 


Assumed  value. 
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Table  10 

Default  parameter  values  for  the  heat  equation 


Symbol 

Quantity 

Size 

k 

Catalyst  layer  thermal  conductivity  [50] 

0.67  Wm-1  K_1 

km 

Membrane  thermal  conductivity  [50] 

0.67  Wm-1  K_1 

kc 

GDL  thermal  conductivity  [50] 

1.67  Wm-1  K-1 

PiQ 

H20  thermal  capacitance 

4.187  xl06Jm-3K-' 

Psci 

Gas  phase  thermal  capacitance/air 

103Jm-3  K-1 

Pm  On 

Membrane  thermal  capacitance3 

2.18  x  106Jm-3K-1 

PcarbQarb 

Carbon  phase  thermal  capacitance 

1.61  x  106Jm“3  K-1 

-Asc 

Entropy  associated  with  ORR  [51  ] 

163.7  J  mol”1  K”1 

-Asa 

Entropy  associated  with  HOR  [51 J 

0J  mol-1  K_1 

3  Estimated. 


values  for  r]b/r^,  r2b/r2j  and  r5f  were  assigned  arbitrarily  and  the 
values  for  ry,  r2j  and  r7f  were  fitted  to  the  provided  a  qualitative 
match  to  the  values  and  trends  in  [1,16].  As  confirmation  of  the 
values,  the  fitted  adsorption  rate  constant  r]  f  was  compared  to  that 
approximated  from  experimental  data  at  different  temperatures  by 
Mohtadi  et  al.  in  [1],  Their  predicted  value  is  of  the  same  order  of 
magnitude,  0(10-2),  as  that  shown  in  Table  6. 

In  the  next  section  use  will  be  made  of  the  spatially  averaged 
reaction  rates  and  surface  coverages.  These  are  defined  as  the  inte¬ 
grals  of  the  reaction  rates  or  surface  coverages  in  space  over  the 


Fig.  2.  The  effect  of  H2S  concentration  and  current  density  on  the  extent  of  poison¬ 
ing  when  operating  in  galvanostatic  mode.  See  Tables  4-10  for  values  of  the  other 
parameters. 


anode  catalyst  layer,  normalised  with  respect  to  L,  the  length  of  the 
catalyst  layer: 


«Ji>=T/  q.'dy 
J  ACL 

<0S>  =  r  f  6s  dy 
J  ACL 


(#H> 


-i/. 


0Hdy. 


The  terminal  value  of  a  quantity  i /r(t)  will  be  defined  as  the  value 
of  iff  at  steady  state,  i.e.: 


“Terminal  value  ofi /r(t)”  =  i/f(oc). 


(41) 


3.  Results  and  discussion 


3.1.  Basic  features  and  validation 

Fig.  2  shows  simulation  results  at  different  current  densities,  at 
a  fixed  H2S  concentration  of  2  ppm,  and  different  H2S  concentra¬ 
tions,  at  a  fixed  current  density  of  0.5  A  cm-2.  The  simulations  are  of 


Fig.  3.  The  effect  of  H2  S  concentration  (top  figure)  and  cell  voltage  (bottom)  on  the 
extent  of  performance  loss  when  operating  in  potentiostatic  mode.  In  the  top  figure 
the  cell  voltage  is  fixed  at  0.5  V  and  in  the  bottom  figure  the  H2S  concentration  is 
fixed  at  2  ppm.  The  temperature  is  60 0  C  in  both  figures.  See  Tables  4-10  for  values 
of  the  other  parameters. 
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operation  in  the  galvanostatic  mode.  The  hydrogen  source  in  these 
simulations  is  40%  H2  in  N2.  The  other  parameter  values  are  given 
in  Tables  4-10  and  are  used  throughout,  unless  otherwise  stated. 
These  tables  also  contain  the  default  values  of  current  density,  volt¬ 
age  and  H2S  concentration.  In  each  simulation  the  initial  conditions 
were  generated  from  a  gradual  increase  in  the  current  over  a  period 
of  500  s  from  OAcm-2,  and  Oppm  H2S,  to  the  stated  current  den- 
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Fig.  4.  The  evolutions  of  the  average  values  of  the  hydrogen  oxidation  rate,  q 4,  and 
the  surface  coverages  of  sulfur  and  hydrogen,  a(q4),  {0s)  and  (0H),  respectively,  at 
different  cell  voltage  values  with  an  H2S  concentration  of  2  ppm  and  at  a  tempera¬ 
ture  of  60  °C.  See  Tables  4-10  for  values  of  the  other  parameter.  Table  11  shows  the 
terminal  values  of  ( 0S )  and  (0H). 


Table  11 

Terminal  values  of  the  spatially  averaged  surface  coverages  of  hydrogen  and  sulfur, 
(0H)  and  (^respectively,  corresponding  to  the  calculations  in  Figs.  3  and  4 

Cell  voltage  (V)  Terminal  (0s)  Terminal  (0h) 

0.3  0.9826  0.01085 

0.5  0.9862  0.00987 

0.7  0.9895  0.00975 

0.9  0.9830  0.01082 


The  H2S  concentration  is  2  ppm  and  the  channel  temperatures  are  60  °C.  See 
Tables  4-10  for  values  of  the  other  parameters. 

sity  value,  followed  by  2  h  at  this  current  density.  A  steady  state  was 
reached  in  each  case. 

In  both  cases,  the  results  in  Fig.  2  results  capture  precisely  the 
trends  observed  in  the  experiments  of  Shi  et  al.  [16],  As  the  current 
density  is  increased  at  fixed  H2  S  concentration  the  extent  of  poison¬ 
ing  (decrease  in  cell  voltage)  increases.  As  the  H2S  concentration  is 
increased  at  a  fixed  current  density  the  degree  of  poisoning  again 
increases. 

Fig.  3  shows  simulation  results  in  the  potentiostatic  mode  at 
different  values  of  H2S  concentration  and  cell  voltage,  again  with 
40%  H2  in  N2.  In  these  simulations  the  initial  conditions  were  gen¬ 
erated  from  a  gradual  increase  in  the  cell  voltage  over  a  period  of 
500  s  from  the  open  circuit  potential,  and  0  ppm  FI2S,  to  the  stated 


Fig.  5.  The  evolutions  of  the  profiles  of  sulfur  and  hydrogen  surface  coverages  during 
poisoning,  corresponding  to  the  calculation  at  a  cell  voltage  of  0.5  V  with  2  ppm  H2S 
in  Fig.  3.  See  Tables  4-10  for  values  of  the  other  parameters. 
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cell  voltage  value,  followed  by  operation  at  this  cell  voltage  until  a 
steady  state  was  reached. 

As  the  H2S  concentration  is  increased  at  a  fixed  cell  voltage  of 
0.5  V,  the  drop  in  current  density  as  a  result  of  poisoning  increases. 
The  profiles  are  a  good  fit  to  the  transient  data  in  [1,18]  showing 
the  decay  of  the  anode  current  density.  The  trend  with  respect 
to  variations  in  cell  voltage  at  a  fixed  H2S  concentration  indicates 
that  the  degree  of  poisoning  initially  increases  as  the  cell  voltage  is 
decreased,  but  as  the  cell  voltage  is  lowered  further,  performance 
improves;  that  is  seen  by  comparing  the  current  density  changes 
from  the  initial  time  to  steady  state,  A i,  shown  in  Fig.  3.  This  result 
is  caused  by  the  trend  in  the  rate  of  hydrogen  oxidation  from  the 
platinum  surfaces,  reaction  (R4),  as  the  cell  voltage  is  lowered.  Fig.  4 
shows  the  evolutions  of  the  spatially  averaged  hydrogen  oxidation 
rate,  a((/4>,  the  spatially  averaged  surface  coverage  of  sulfur,  (6s), 
and  the  spatially  averaged  surface  coverage  of  hydrogen,  (0H>-  In 
these  figures  the  cell  voltage  is  varied  from  0.9  to  0.3  V  with  2  ppm 
H2S  at  an  operating  temperature  of  60  °C,  as  in  Fig.  3. 

The  variation  in  current  density  in  Fig.  3  as  the  cell  voltage  is 
lowered  mirrors  the  variation  in  a(q4),  clearly  demonstrating  that 
the  oxidation  of  hydrogen  controls  the  extent  of  degradation.  The 
terminal  values  of  (6s)  and  (%),  defined  in  Eq.  (41),  are  shown  in 
Table  11.  The  terminal  value  of  (#h)  decreases  as  the  cell  voltage 
is  increased  from  0.3  to  0.5  V  and  again  from  0.5  to  0.7  V.  On  the 
other  hand,  the  terminal  value  of  (#h>  increases  as  the  cell  voltage 
is  increased  from  0.7  to  0.9  V.  Correspondingly,  the  terminal  value 
of  (6s)  increases  from  0.3  to  0.7  V  and  decreases  from  0.7  to  0.9  V. 
This  would  suggest  that  as  the  cell  voltage  is  lowered  the  extent 
of  degradation,  measured  by  the  drop  in  current  density,  reaches  a 
maximum,  which  is  indeed  the  case  in  Fig.  3. 

Fig.  5  shows  the  evolution  of  the  surface  coverage  profiles  during 
the  calculation  at  0.5  V  with  2  ppm  H2S.  At  the  initial  time,  when 
the  concentration  of  H2S  is  increased  from  0  to  2  ppm,  the  coverage 
of  platinum  sites  by  hydrogen  is  almost  complete.  Within  2.78  h 
roughly  40%  of  the  surface  area  is  covered  by  sulfur,  with  a  corre¬ 
sponding  decrease  in  adsorbed  hydrogen.  Thereafter,  the  adsorp¬ 
tion  rate  of  sulfur  is  increasingly  slower  but  at  the  terminal  current 
density  the  surface  coverage  of  sulfur  is  close  to  a  monolayer:  6S  1 . 

3.2.  Effect  of  channel  temperature  and  water  activity 
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Fig.  7.  The  evolution  of  the  spatially  averaged  surface  coverages  of  sulfur  and 
hydrogen,  {0s)  and  (On),  respectively,  corresponding  to  the  calculations  at  channel 
temperatures  of  60  and  80  °  C  shown  in  Fig.  6.  See  Tables  4-10  for  values  of  the  other 
parameters. 


The  effect  of  temperature  on  the  extent  of  poisoning  is  depicted 
in  Fig.  6,  showing  simulation  results  at  different  channel  tempera- 
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Fig.  6.  The  effect  of  channel  temperature  on  H2  S  poisoning:  for  an  H2  S  concentration 
of  2  ppm  and  a  cell  voltage  of  0.5  V.  See  Tables  4-10  for  values  of  the  other  parameters. 


tures  for  an  H2S  concentration  of  2  ppm  and  a  cell  voltage  of  0.5  V.  It 
is  immediately  noticeable  that  the  initial  current  density  increases 
as  the  channel  temperature  is  increased.  There  are  several  effects 
associated  with  a  temperature  increase  at  fixed  values  of  other 
quantities  in  the  absence  of  H2S:  a  decrease  in  the  channel  con¬ 
centration  of  reactants:  a  decrease  in  the  reaction  rates  for  oxygen 
reduction  and  hydrogen  oxidation  through  the  Arrhenius  depen¬ 
dence;  a  decrease  in  the  rate  of  condensation  through  an  increase 
in  the  saturation  vapour  pressure  and  a  reduction  in  the  channel 
vapour  concentration;  and  an  increase  in  the  membrane  conductiv¬ 
ity.  The  latter  effect  dominates,  leading  to  the  higher  initial  current 
density. 

As  the  channel  temperature  is  increased  the  drop  in  current 
density  from  its  initial  to  terminal  value  decreases,  leading  to  a 
reduced  degree  of  poisoning.  This  can  bee  seen  from  the  current 
density  changes,  A i,  shown  in  Fig.  6  at  each  temperature.  Fig.  7 
shows  the  evolution  of  the  spatially  averaged  surface  coverages  of 
sulfur  and  hydrogen  defined  in  Eq.  (40),  corresponding  to  the  cal¬ 
culations  at  60  and  80 0  C  in  Fig.  6.  These  plots  demonstrate  that 
the  reduced  degree  of  poisoning  is  the  result  of  a  greater  terminal 
coverage  of  hydrogen  and  a  lower  terminal  coverage  of  sulfur  at  the 
higher  temperature. 
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However,  Fig.  7  also  shows  that  the  coverage  of  platinum  sites 
by  sulfur  is  more  rapid  as  the  temperature  is  increased,  so  that  0$ 
is  higher  and  9»  is  lower  up  to  approximately  t  =  13  h,  at  which 
point  the  curves  intersect.  Thus,  although  poisoning  is  less  severe 
at  the  higher  temperature  on  a  long  time  scale,  the  rate  of  poisoning 
is  faster,  i.e.,  the  current  density  drops  more  rapidly.  Indeed  the 
two  curves  at  60  and  80°  C  in  Fig.  6  intersect  at  approximately 
t  =  7h. 

The  above  result  agrees  with  the  experimental  data  of  Mohtadi 
et  al.  [1  ],  in  which  the  authors  estimate  the  surface  coverage  of  sul¬ 
fur  at  different  temperatures  by  dividing  the  oxidation  charge  by  the 
maximum  total  charge  obtained  under  sulfur  saturation  conditions. 
More  precisely,  the  shape  of  the  profiles  of  (9s)  and  the  trend  with 
respect  to  temperature  in  Fig.  7,  before  a  steady  state  is  established, 
match  the  shape  and  trend  in  Fig.  4  of  [1].  The  timescales  are  also 
very  similar.  Mohtadi  et  al.  do  however  estimate  that  the  surface 
coverage  of  sulfur  at  long  times  increases  with  temperature,  which 
seems  not  to  support  the  crossing  of  the  curves  in  Fig.  7.  The  authors 
also  state  that  the  half-cell  current  decrease  at  90  °  C  was  lower 
than  the  current  decrease  corresponding  to  the  same  sulfur  cover¬ 
age  for  full  cells  at  70  and  50  °C.  Therefore  they  chose  not  to  base 
the  coverage  calculations  on  the  current  decrease.  They  hypothesise 
that  sulfur  crosses  over  to  the  cathode,  thus  affecting  the  oxygen 
reduction  reaction.  Furthermore,  the  full-cell  curves  at  50  and  70  ° 
C  in  Fig.  7  of  their  paper  were  found  to  cross  at  intermediate  times 
before  a  steady  state  was  established,  in  a  manner  similar  to  that 
in  Fig.  6.  The  results  are  therefore  inconclusive.  Certain  elements 
of  the  results  displayed  in  Figs.  6  and  7  can  be  validated,  but  the 
behaviour  seen  as  the  steady  states  are  approached  cannot  be  con¬ 
firmed  until  further  data  is  available.  On  the  other  hand  several 
of  the  results  in  [1]  seem  to  support  it.  The  following  hypothesis, 
to  be  discussed  below,  is  made:  the  nonlinear  behaviour  found  in 


[1]  is  due  to  the  competition  between  adsorption  and  oxidation  of 
sulfur. 

To  investigate  the  behaviour  described  above  the  reader  is 
referred  to  Fig.  8,  which  shows  the  corresponding  evolutions  at 
both  temperatures  of  the  spatially  averaged  rates  of  the  surface 
reactions  (Rl),  (R3),  (R4)  and  (R7),  defined  in  Eq.  (40):  a(q\),  a(q3), 
a(q 4>  and  0(137).  These  are  sulfur  adsorption,  hydrogen  adsorption, 
hydrogen  oxidation  and  sulfur  oxidation  respectively.  In  both  cases 
the  average  rate  of  sulfur  adsorption,  a(q j),  rises  rapidly  as  the  fuel 
source  is  changed  from  clean  to  contaminated.  The  terminal  value 
of  a(q  i )  increases  as  the  temperature  is  increased,  and  a(q,}  relaxes 
to  a  steady-state  value  in  a  shorter  period  of  time.  At  fixed  9$  and 
0pt,  the  forward  part  of  the  reaction  rate  q\  in  Eq.  (29)  increases 
and  the  backward  part  decreases  as  temperature  is  increased.  The 
result  is  a  higher  value  of  9s  at  80 0  C  during  the  early  stages  of 
its  evolution  (before  steady  state).  This  occurs  in  competition  with 
the  increased  adsorption  rate  of  hydrogen  shown  in  Fig.  8.  As  9s 
increases  the  sulfur  adsorption  rate  decreases  (a  consequence  of 
the  Frumkin  kinetics),  and  thus  a  steady  state  is  reached  sooner.  On 
the  other  hand,  the  average  oxidation  rate  of  sulfur,  a(q7),  is  greater 
at  the  higher  temperature.  The  relative  increase  in  the  terminal 
value  of  a(q7),  that  is 

terminal  value  at  80  °C  —  terminal  value  at  60  °C 
terminal  value  at  60  °C 

is  equal  to  1.9167,  while  the  relative  increase  in  the  terminal  value 
of  a(q\ )  is  1.727,  which  further  suggests  that  the  terminal  value  of 
9s  should  be  lower  at  80  °C.  This  is  indeed  the  case  in  Fig.  7.  In 
other  words,  the  increase  in  the  adsorption  rate  of  sulfur  at  80 0 
C  is  dominated  by  the  simultaneous  increase  in  its  oxidation  rate. 
The  equivalent  relative  increases  in  a(q3)  and  a(q4)  are  almost  iden¬ 
tical  (approximately  0.62)  implying  that  the  increased  adsorption 


Fig.  8.  The  evolution  of  the  spatially  averaged  rates  of  the  surface  reactions  (Rl ),  (R3),  (R4)  and  (R7),  a(<ji ),  a{q3),  a{q4)  and  a(q7),  respectively,  for  the  1 
shown  in  Fig.  6.  See  Tables  4-10  for  values  of  the  other  parameters. 
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Fig.  9.  The  effect  of  simultaneous  changes  in  the  channel  temperature  and  cell  volt¬ 
age  onH2S  poisoning:  for  anH2S  concentration  of  2  ppm.  See  Tables  4-10  for  values 
of  the  other  parameters. 


of  atomic  hydrogen  is  dominated  by  its  increased  oxidation.  The 
higher  terminal  value  of  0H  at  80°  C  is  therefore  a  direct  con¬ 
sequence  of  the  lower  terminal  value  of  6s.  Note  finally  that  the 
reaction  rates  for  sulfur  are  much  smaller  than  those  for  hydrogen, 
which  is  due  to  the  very  small  H2S  concentrations. 

These  results  pertain  to  a  cell  voltage  of  0.5  V.  To  see  the  effect 
of  changes  in  the  cell  voltage  as  the  channel  temperature  is  var¬ 
ied,  the  reader  is  referred  to  Fig.  9.  As  the  cell  voltage  is  increased, 
the  difference  in  the  terminal  current  density  between  60  and  80 0 
C  decreases.  For  example,  at  £  =  0.3  V  the  difference  is  approxi- 
mately0.37  Acm~2andat£  =  0.7Vitis  approximately  0.14  A  cm-2. 
Moreover,  there  is  a  qualitative  difference  between  the  two  sets  of 
curves  in  Fig.  9;  at  £  =  0.3  V  performance  is  uniformly  better  at 
higher  temperature,  whereas  at  £  =  0.7  the  performance  at  lower 
temperature  is  superior  except  towards  the  beginning  of  the  calcu¬ 
lation  and  towards  the  end  after  the  terminal  current  density  has 
been  reached  at  the  higher  temperature;  performance  at  60°  C  is 
better  than  performance  at  70 0  C  between  3 .39  <  t  <  16.26  h.  For 
both  £  =  0.3  V  and  £  =  0.7  V,  Fig.  10  shows  the  evolutions  of  a(q\) 
and  0(97)  at  80°  C  and  Fig.  11  shows  the  evolutions  of  a(<ji >  and 
a(q7)  at  60 °C.  At  60  °C,  there  is  a  visible  delay  in  the  times  taken 
for  both  a(q, )  and  a(q7 )  to  relax  to  their  steady-state  values.  In  con¬ 
trast,  the  differences  in  these  times  at  80 0  C  are  slight.  Thus,  while 


at  80 0  C  the  rate  of  poisoning  changes  little  when  the  cell  voltage  is 
decreased  from  0.7  to  0.3  V,  at  60 0  C  it  visibly  increases.  This,  allied 
with  the  increasing  differences  in  initial  current  densities  as  the 
cell  voltage  is  increased,  explains  why  the  curves  in  Fig.  9  intersect 
at  0.7  V  but  not  at  0.3  V. 

One  of  the  main  performance  control  mechanisms  in  PEM  fuel 
cells  is  the  channel  water  activity,  which  must  be  chosen  such  that 
the  membrane  remains  well  hydrated  for  sufficiently  high  protonic 
conductivity  and  such  that  flooding  does  not  occur  in  the  cathode 
catalyst  layer,  restricting  reactant  access  to  the  catalyst  sites.  Fig.  12 
shows  the  effect  of  variations  in  the  water  activity  (the  same  in  both 
channels)  at  a  cell  voltage  of  0.5  V  and  with  2  ppm  H2 S  in  40%  H2 /N2 . 
All  other  parameters  are  fixed  as  in  Tables  4-10.  The  plots  reveal  an 
improvement  in  performance  in  going  from  aw  =  0.7  (equilibrium 
relative  humidity  of  70%)  to  aw  =  1  in  both  channels.  However,  for 
water  activities  below  0.7,  very  little  change  in  the  performance  is 

It  is  instructive  to  examine  the  water  saturation  profiles  for  the 
two  cases  ofaw  =  1  and  aw  =  0.7;  these  are  shown  in  Fig.  ^.Imme¬ 
diately  noticeable  is  that  the  saturation  levels  in  the  case  aw  =  0.7 
are  lower,  clearly  as  a  result  of  the  reduced  rates  of  condensation. 
The  levels  in  the  anode  are  practically  zero  and  are  not  therefore 
visible  in  these  plots.  The  second  feature  to  notice  is  that  the  sat¬ 
uration  levels  fall  markedly  during  the  poisoning  process.  In  the 
cathode  the  lower  current  density  due  to  poisoning  reduces  the 


Fig.  10.  The  evolutions  of  the  spatially  averaged  rates  of  the  surface  reactions  (R1 ) 
and  (R7),  alq, )  and  a(q7),  respectively,  for  the  case  80°  C  shown  in  Fig.  9.  See 
Tables  4-10  for  values  of  the  other  parameters. 
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Fig.  11.  The  evolutions  of  the  spatially  averaged  rates  of  the  surface  reactions  (R1 ) 
and  (R7),  a(qj)  and  a(q7),  respectively,  for  the  case  60°  C  and  2  ppm  H2S.  See 
Tables  4-10  for  values  of  the  other  parameters. 


Fig.  12.  The  effect  of  channel  water  activity  on  H2S  poisoning  for  an  H2S  concen¬ 
tration  of  2  ppm  and  a  cell  voltage  of  0.5  V.  See  Tables  4-10  for  values  of  the  other 
parameters. 


activities  of  1  and  0.7  shown  in  Fig.  12.  See  Tables  4-10  for  values  of  the  other 
parameters.  The  left-hand  boundary,  y  =  0,  corresponds  to  the  cathode  channel/gas 
diffusion  layer  interface  and  the  right-hand  boundary,  y  =  500  p,m,  to  the  anode 
channel/gas  diffusion  layer  interface.  The  region  225  p,m<y<275  p,m  corresponds 
to  the  membrane. 

rate  of  water  production  and  in  the  anode  water  is  consumed  in  the 
oxidation  reaction  (R7).  The  lower  current  density  will  also  reduce 
back  diffusion  caused  by  the  drag  of  water  molecules  attached 
to  protons,  from  the  anode  to  cathode,  as  is  seen  from  Eq.  (6).  A 
main  consequence  of  the  reduced  saturation  levels  is  a  decrease 
in  the  membrane  conductivity,  which  depends  on  the  water  con¬ 
tent  in  the  membrane  according  to  Eq.  (14).  The  water  content  in 
turn  is  a  function  of  the  local  saturation  and  water  vapour  levels. 
Thus,  the  reduction  in  current  density  due  to  sulfur  coverage  of  the 
platinum  sites  is  compounded  by  a  reduced  membrane  conductiv¬ 
ity. 

4.  Conclusions 

A  modelling  framework  for  predicting  and  studying  the  poi¬ 
soning  effect  of  H2S  on  the  anode  catalyst  layer  of  a  PEM  fuel  cell 
has  been  developed.  In  contrast  to  the  model  in  [16],  the  present 
model  explicitly  include  mass,  energy  and  momentum  conserva¬ 
tion,  together  with  the  fundamental  modes  of  transport  and  a  more 
detailed  kinetic  mechanism.  The  model  can  simulate  both  galvano- 
static  and  potentiostatic  operation.  Comparison  with  data  available 
in  the  open  literature  has  shown  that  the  trends  are  well  captured 
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with  respect  to  variations  in  the  H2S  concentration,  cell  voltage  and 
current. 

The  model  predicts  that  temperature  and  water  activity  vari¬ 
ations  have  a  complicated  effect  on  the  poisoning  process.  The 
kinetic  mechanism  in  the  anode  is  intimately  linked  with  varia¬ 
tions  in  these  quantities,  which  also  yield  a  wider  influence  on 
performance,  through  the  form  of  the  reaction  rates.  For  exam¬ 
ple,  oxidation  of  sulfur  from  the  platinum  surfaces  in  the  anode 
consumes  water,  and  is  therefore  dependent  on  the  channel  water 
activity.  It  has  been  demonstrated  that  the  anode  water  levels  can 
decrease  quite  significantly  as  a  consequence  of  this  reaction.  The 
decreased  water  levels  in  the  anode  will  reduce  the  membrane 
conductivity,  further  decreasing  the  current  density.  The  reduced 
current  density  as  a  result  of  poisoning  in  turn  reduces  the  water 
production  rate  in  the  cathode  and  restricts  back  diffusion  of  water 
via  proton  migration.  The  relative  strengths  of  these  effects  can  be 
investigated  with  the  model. 

Temperature  increases  were  shown  in  general  to  lessen  the 
degree  of  poisoning,  although  the  simulation  results  also  suggest 
that  the  behaviour  of  the  system  over  a  relatively  short  timescale 
(before  any  steady  state  is  reached)  is  not  straightforward  to  pre¬ 
dict  and  certain  features  may  be  masked  by  the  steady-state  results. 
For  these  relatively  short  timescales  the  results  of  this  study  agree 
with  those  of  Mohtadi  et  al.  [  1  ],  but  at  steady  state  the  comparison  is 
inconclusive.  It  has  been  hypothesised  that  the  discrepancies  found 
by  Mohtadi  et  al.  could  be  due  to  the  relative  changes  in  the  rates 
of  sulfur  adsorption  and  oxidation  as  the  operating  temperature  is 
varied.  The  adsorption  rate  almost  certainly  increases  as  the  tem¬ 
perature  is  raised,  but  one  can  predict  that  there  is  a  corresponding 
increase  in  the  oxidation  rate,  which  over  long  times  could  lead  to 
better  performance  at  higher  temperature.  More  data  is  needed  to 
verify  these  claims. 

In  this  work  several  assumptions  have  been  employed  in  relation 
to  the  kinetic  model.  The  missing  or  estimated  electrochemical  con¬ 
stants  can  be  estimated  from  experiment  and  input  to  the  model, 
improving  its  accuracy.  Future  work  will  focus  on  obtaining  these 
parameters  and  validating  the  model  over  a  much  broader  range  of 
operating  conditions. 
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